Analysis of Standard Morbidity Ratio on Recorded Dry Season Malaria Incidences in Kasungu
Clinton Nkolokosa
Last edited 03 June 2021
This geospatial statistical model uses routinely collected malaria case data, population data and remotely sensed data, such as open and vegetated water bodies, to estimate population living around open water bodies, and ultimately quantify the association between proximity to larval habitat and malaria risk in health facility catchment areas in Kasungu. Buffers around waterbodies are created and then combined with population data in raster format to estimate the total population that is within the buffer. Subsequently, the observed malaria cases are modeled using Poisson regression to find out if living within various distances from water bodies is causing variability in malaria risk in Kasungu district. We hypothesize that the risk of being a case in a catchment is dependent on proximity to water bodies. The data used spans from 2015 to 2020 and was derived from digitized DHIS2 malaria records, accessibility mapping, aggregated population geospatial layer and TropWet tool in Google Earth Engine.
Load packages
Loading the R packages that will be used to read in, view, transform and model the malaria cases and spatial datasets.
library(spatialEco)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(plotly)
library(lubridate)
library(knitr)
library(raster)
library(rgdal)
library(rgeos)
library(sf)
library(sp)
library(tmap)
library(spdep)
library(maptools)
library(gridExtra)
library(grid)
library(exactextractr)
library(DataExplorer)
library(mapview)
`%>%` <- magrittr::`%>%`Tell R where the data is
here::here()## [1] "C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021"Load datasets
The total dry season malaria cases recorded at health-care facilities in Kasungu from 2017 to 2019 are contained in the KasunguData.csv sourced from https://dhis2.health.gov.mw/. The kasungu_facility_catchments_2004.shp shapefile also contains the population and health information within each health-facility catchment area in Kasungu district. The aggregated population raster layers for Malawi e.g.,ku_pop_2017_1km_aggregated.tif were downloaded from the Open Spatial and Demographic and Data Research website: https://www.worldpop.org/geodata/country?iso3=MWI. These layers estimate total number of people per grid-cell. The units are number of people per pixel with country totals adjusted to match the corresponding official United Nations population estimates. The datasets were downloaded in Geotiff at a resolution of 1km and are projected in Geographic Coordinate System, WGS84. The kasungu_water.shpand water_bodies layers contain open and vegetated waterbodies polygons, detected using the Tropical Wetland Unmixing Tool (TropWet). TropWet is a Google Earth Engine hosted toolbox that uses the Landsat archive to map tropical wetlands and can be accessed through: https://www.aber.ac.uk/en/dges/research/earth-observation-laboratory/research/tropwet/
# 2017, 2018 and 2019 dry season malaria cases
dry_season_malaria_2015_2019 <- read.csv(here::here("data", "dry_season_malaria 2015-2019.csv"))
# Export 2017, 2018 and 2019 dry season malaria cases
write.csv(dry_season_malaria_2015_2019, file = "data/dry_season_malaria_2017_2019.csv")
write.table(dry_season_malaria_2015_2019,
file = "dry_season_malaria_2017_2019.txt",
sep = ",",
quote = FALSE,
row.names = FALSE)
# 2020 dry season malaria cases
ku_malaria_2020 <- read.csv(here::here("data", "dry_season_malaria_2020.csv"))
# Merge 2015 to 2019 dry season malaria case data with 2020 data
dry_season_malaria_2017_2020 <- cbind.data.frame(dry_season_malaria_2015_2019, ku_malaria_2020)
dry_season_malaria_2017_2020 <- dry_season_malaria_2017_2020[,c("FID", "Names", "dr_2017",
"dr_2018", "dr_2019", "dr_2020",
"LONGITU", "LATITUD")]
# Export 'dry season malaria 2017-2020' as csv
write.csv(dry_season_malaria_2017_2020, file = "data/dry_season_malaria_2017_2019.csv")
# Kasungu district boundary shapefile
kasungu_district <- st_read(here::here("data", "kasungu_district.shp"))## Reading layer `kasungu_district' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_district.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 491272.7 ymin: 8494349 xmax: 609044.2 ymax: 8632164
## Projected CRS: WGS 84 / UTM zone 36S# Kasungu health facility catchments enerated from accessibility mapping
malire_new <- sf::st_read(here::here("data", "Kasungu_new_health_fac_catchment_clip.shp")) ## Reading layer `Kasungu_new_health_fac_catchment_clip' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\Kasungu_new_health_fac_catchment_clip.shp' using driver `ESRI Shapefile'
## Simple feature collection with 29 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 491878.5 ymin: 8494645 xmax: 609044.2 ymax: 8631908
## Projected CRS: WGS 84 / UTM zone 36S# Kasungu population raster layer
kasungu_population_2017 <- raster(here::here("data", "ku_pop_2017_1km_aggregated.tif"))
kasungu_population_2018 <- raster(here::here("data", "ku_pop_2018_1km_aggregated.tif"))
kasungu_population_2019 <- raster(here::here("data", "ku_pop_2019_1km_aggregated.tif"))
kasungu_population_2020 <- raster(here::here("data", "ku_pop_2020_1km_aggregated.tif"))
# Open waterbodies polygons
water_2017 <- st_read(here::here("data", "water_bodies_2017.shp"))## Reading layer `water_bodies_2017' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\water_bodies_2017.shp' using driver `ESRI Shapefile'
## Simple feature collection with 168 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 514497 ymin: 8495941 xmax: 603149.8 ymax: 8620169
## Projected CRS: WGS 84 / UTM zone 36Swater_2018 <- st_read(here::here("data", "kasungu_2018_water.shp"))## Reading layer `kasungu_2018_water' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_2018_water.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1105 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 496807.6 ymin: 8494693 xmax: 607913.8 ymax: 8607747
## Projected CRS: WGS 84 / UTM zone 36Swater_2019 <- st_read(here::here("data", "kasungu_2019_water.shp"))## Reading layer `kasungu_2019_water' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_2019_water.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1941 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 494197.2 ymin: 8494693 xmax: 607913.8 ymax: 8617573
## Projected CRS: WGS 84 / UTM zone 36Swater_2020 <- st_read(here::here("data", "water_bodies_2020.shp"))## Reading layer `water_bodies_2020' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\water_bodies_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 266 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 508985.6 ymin: 8495793 xmax: 585761.1 ymax: 8620169
## Projected CRS: WGS 84 / UTM zone 36S# Add a field ID to water bodies polygons
water_2017$ID <- 1:nrow(water_2017)
water_2018$ID <- 1:nrow(water_2018)
water_2019$ID <- 1:nrow(water_2019)
water_2020$ID <- 1:nrow(water_2020)View the dry season malaria case data
We observe that Kasungu district has 30 health facilities classified as dispensary, health centre, district hospital and rural hospital and the highest malaria cases were recorded at Kasungu District Hospital.
dry_season_malaria_2017_2020 %>%
glimpse() %>%
summary()## Rows: 36
## Columns: 8
## $ FID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,~
## $ Names <chr> "Anchor Farm", "Bua Health Centre", "Chamama Health Facility",~
## $ dr_2017 <int> 2966, 3489, 1878, NA, 3601, 2292, 5801, 2192, 2745, NA, NA, NA~
## $ dr_2018 <int> 3142, 1804, 1750, NA, 4027, 2116, 4502, 2394, 4138, NA, NA, NA~
## $ dr_2019 <int> 1978, 1740, 1508, NA, 1686, 1770, 5470, 2553, 2200, NA, NA, NA~
## $ dr_2020 <int> 2837, 6194, 1490, 542, 1386, 2411, 9628, 4196, 3957, 147, 3003~
## $ LONGITU <dbl> 33.39000, 33.53859, 33.77686, NA, 33.68528, 33.79701, 33.30816~
## $ LATITUD <dbl> -13.41000, -13.29403, -12.92431, NA, -13.10757, -12.97452, -12~## FID Names dr_2017 dr_2018
## Min. : 1.00 Length:36 Min. : 427 Min. : 749
## 1st Qu.: 9.75 Class :character 1st Qu.: 2196 1st Qu.: 2424
## Median :18.50 Mode :character Median : 3132 Median : 3357
## Mean :18.50 Mean : 3586 Mean : 3909
## 3rd Qu.:27.25 3rd Qu.: 3993 3rd Qu.: 4684
## Max. :36.00 Max. :16289 Max. :15821
## NA's :6 NA's :6
## dr_2019 dr_2020 LONGITU LATITUD
## Min. : 533 Min. : 0 Min. :33.18 Min. :-13.57
## 1st Qu.: 1748 1st Qu.: 2286 1st Qu.:33.38 1st Qu.:-13.25
## Median : 2556 Median : 3824 Median :33.50 Median :-12.98
## Mean : 2880 Mean : 4657 Mean :33.52 Mean :-12.99
## 3rd Qu.: 3284 3rd Qu.: 6212 3rd Qu.:33.68 3rd Qu.:-12.79
## Max. :10721 Max. :24424 Max. :33.87 Max. :-12.42
## NA's :6 NA's :6 NA's :6dry_season_malaria_2017_2020 %>%
plotly::plot_ly(y = ~Names,
x = ~dr_2017,
type = "bar",
orientation = 'h',
name = "2017") %>%
plotly::add_trace(x = ~ dr_2018,
name = "2018") %>%
plotly::add_trace(x = ~ dr_2019,
name = "2019") %>%
plotly::add_trace(x = ~ dr_2020,
name = "2020") %>%
plotly::layout(#barmode = "stack",
xaxis = list(title = "Total malaria cases"),
yaxis = list(title = ""),
hovermode = "compare",
margin = list(b = 10,
t = 10,
pad = 2))Fig.1 The total malaria cases recorded at each health-care facility in Kasungu district
View location of Kasungu health-care facilities within health facility catchment area
Heath facility catchment area is the area from which a health facility attracts patients. The new health facility catchments polygon was generated from generic accessibility mapping script adapted from https://malariaatlas.org/wp-content/uploads/accessibility/R_generic_accessibilty_mapping_script.r The script requires two user supplied datasets: the 2015 friction surface, which is available here: http://www.map.ox.ac.uk/accessibility_to_cities/, and a user-supplied .csv of points dry_season_malaria_2017_2020. The accumulated cost algorithm accCost and r.Cost algorithm in GRASS GIS were run to make the final output map of new health facility catchment boundaries.
# Using the complete.cases() function to select health centres with complete longitude and latitude coordinates
# Aggregating Kasalika Health Centre and Kasungu District Hospital, and
# Kaluluma Rural Hospital and Nkhamenya Rural Hospital malaria cases
zipatala_aggregated <- dry_season_malaria_2017_2020[complete.cases(dry_season_malaria_2017_2020),] %>%
dplyr::filter(Names != "Kasalika Health Centre",
Names != "Bua Health Centre",
Names != "Kaluluma Rural Hospital")
zipatala_aggregated$dr_2017[which(zipatala_aggregated$Names == "Kasungu District Hospital")] <- 4528 + 16289
zipatala_aggregated$dr_2018[which(zipatala_aggregated$Names == "Kasungu District Hospital")] <- 4493 +15821
zipatala_aggregated$dr_2019[which(zipatala_aggregated$Names == "Kasungu District Hospital")] <- 2729 + 10721
zipatala_aggregated$dr_2020[which(zipatala_aggregated$Names == "Kasungu District Hospital")] <- 4368 + 24424
zipatala_aggregated$dr_2017[which(zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- 2887 + 752
zipatala_aggregated$dr_2018[which(zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- 851 + 3689
zipatala_aggregated$dr_2019[which(zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- 533 + 4004
zipatala_aggregated$dr_2020[which(zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- 3587 + 5929
zipatala_aggregated$dr_2017[which(zipatala_aggregated$Names == "Mziza Health Centre")] <- 3863 + 3489
zipatala_aggregated$dr_2018[which(zipatala_aggregated$Names == "Mziza Health Centre")] <- 2815 + 1804
zipatala_aggregated$dr_2019[which(zipatala_aggregated$Names == "Mziza Health Centre")] <- 2439 + 1740
zipatala_aggregated$dr_2020[which(zipatala_aggregated$Names == "Mziza Health Centre")] <- 6194 + 2397
health_facility_aggr_sf <- sf::st_as_sf(zipatala_aggregated,
coords = c("LONGITU", "LATITUD"),
crs = 4326, agr = "constant")
tm_shape(malire_new)+
tm_polygons()+
tm_shape(health_facility_aggr_sf)+
tm_dots(size = .3,
col = "blue",
alpha = 0.5)+
tm_text("Names",
size = .3,
just = "top",
col = "black",
remove.overlap = TRUE)+
tm_layout(frame = FALSE,
title = "New Kasungu health facility \n catchment boundaries",
title.size = .8,
title.position = c("left", "top"))+
tm_compass(position=c("right", "top"))+
tm_scale_bar(breaks = c(0, 10, 20),
text.size = .5)
Fig 2. Kasungu Health-care Facilities and Catchment Areas
View raster population metadata and estimated population per grid-cell
# Take a glimpse at the WorldPop raster layers
kasungu_population_2017## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2017_1km_aggregated.tif
## names : ku_pop_2017_1km_aggregatedkasungu_population_2018## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2018_1km_aggregated.tif
## names : ku_pop_2018_1km_aggregated
## values : 0, 6253.557 (min, max)kasungu_population_2019## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2019_1km_aggregated.tif
## names : ku_pop_2019_1km_aggregated
## values : 0, 6483.727 (min, max)kasungu_population_2020## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2020_1km_aggregated.tif
## names : ku_pop_2020_1km_aggregated
## values : 0, 7949.033 (min, max)# Function to create a raster population map
create.population.map <- function(population.raster, title){
# raster population map
# arguments:
# population.raster: aggregated population raster layer from WorldPop
# legend.title: legend title
# returns:
# a tmap-element (plots a map)
tm_shape(population.raster)+
tm_raster(palette = "YlOrBr",
title = title,
breaks = c(0,100,200,400,600,800,1000,2000,4000,6000,8000))+
tm_layout(legend.position = c("right", "bottom"),
frame = FALSE)+
tm_scale_bar(position = c("left", "bottom"))
}
# Set to static map
tmap_mode("plot")
estimated_pop_2017 <- create.population.map(kasungu_population_2017, title = "2017 Population")
estimated_pop_2018 <- create.population.map(kasungu_population_2018, title = "2018 Population")
estimated_pop_2019 <- create.population.map(kasungu_population_2019, title = "2019 Population")
estimated_pop_2020 <- create.population.map(kasungu_population_2020, title = "2020 Population")
# Layout the maps
tmap_arrange(estimated_pop_2017, estimated_pop_2018, estimated_pop_2019, estimated_pop_2020, nrow = 2) 
Fig.4 Estimated total number of people per 1km grid-cell
Assign dry season malaria cases and population density to new health facility catchments
The WorldPop aggregated population e.g. kasungu_population_2017.tif, and DHIS2 malaria data dry_season_malaria_2017_2020 datasets on population and malaria are assigned to the new health facility catchments.
# Create a function that assigns malaria data from health facilities to their catchments areas --------------------
assign.malaria.data <- function(catchment_boundary, malaria_data){
# function to assign malaria cases from health facilities to their corresponding catchment
# arguments:
# catchment_boundary: sf polygon object of new catchment boundaries
# malaria_data: sf point object with a data frame containing the dry season malaria cases
# returns:
# catchments_malaria_sf: sf polygon object with a data frame containing dry season malaria cases
# Convert sf objects to spatial
catchment_shp <- as(catchment_boundary, "Spatial")
malaria_shp <- as(malaria_data, "Spatial")
# Match CRS
malaria_shp <- spTransform(malaria_shp, crs(catchment_shp))
# Overlay aggregated health facility points and extract 2017 - 2020 malaria cases
# Using 'point.in.poly' to return a point spatial object, in this case location of health facilities
# and estimated population instead of sp::over function, which simply returns
# a data frame, with the same no. rows.
# Argument 'sp = TRUE' returns an sp class object, else returns sf class object
# Joining the malaria and population dataset using only 'merge' function can't work due to
# non-unique columns and differences in row numbers
hospitals_in_catchment <- spatialEco::point.in.poly(malaria_shp, catchment_shp, sp = TRUE)
# Add the extracted ID, health facility names and dry season malaria cases to
# the health facility catchments (hfc)
hfc_malaria_shp <- merge(catchment_shp, hospitals_in_catchment,
by.x = "DN", by.y = "FID")
# Convert the shapefile containing malaria data to sf-object
hfc_malaria_sf <- sf::st_as_sf(hfc_malaria_shp)
# Tidy the data by dropping columns not needed
catchment_malaria <- hfc_malaria_sf %>%
dplyr::select(-c(coords.x1, coords.x2))
return(out = catchment_malaria)
}
# Invoking the function ----------------------------------------------------------------------------------
malaria_by_catchment <- assign.malaria.data(malire_new, health_facility_aggr_sf)
# Create a function to extract population from WroldPop raster file and assign ---------------------------
# the values to the new catchments.
extract.pop.values <- function(kasungu_pop_raster, catchments){
# function to extract population from raster file and assign the population to catchments
# arguments:
# kasungu_pop_raster: population raster file clipped to Kasungu district
# catchments: shapefile containing the polygons that we wish to use as boundaries
# returns:
# catchments_malaria_pop_sf: sf polygon object with a data frame containing malaria and population data
# Crop and mask the population raster to exclude Kasungu National Park
pop_raster_clip <- raster::mask(raster::crop(kasungu_pop_raster, extent(catchments)), catchments)
# Extracting zonal statistics from a population raster layer.
# The population raster is a continuous gridded surface layer that has an
# estimated population density value to every square in their grid.
# The population values are then summed and apportioned to the catchment polygons
# catchments_malaria_pop <- catchments %>%
# dplyr::mutate(pop = round(raster::extract(pop_raster_clip, catchments, fun = sum, na.rm = TRUE)))
pop_by_catchment <- round(raster::extract(pop_raster_clip, catchments, fun = sum, na.rm = TRUE))
pop_by_catchment_df <- pop_by_catchment %>%
# apply unlist to the lists to have vectors as the list elements
lapply(unlist) %>%
# convert vectors to data.frames
lapply(as_tibble) %>%
# combine the list of data.frames
bind_rows(., .id = "rowID") %>%
# rename the value variable
dplyr::rename(pop = value)
# Add row ID to column to catchment layer
catchments$rowID <- 1:nrow(catchments)
# Merge catchment areas and population data
pop_by_catchments <- merge(catchments, pop_by_catchment_df, by.x = "rowID", by.y = "rowID")
# Reorder columns by position
# Get column names colnames()
# [1] "rowID" "DN" "FID" "Names" "dr_2017" "dr_2018" "dr_2019" "dr_2020" "pop" "geometry"
# pop_by_catchments <- pop_by_catchments[, c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)]
# Cleaning 'Inf' values
pop_by_catchments %>%
dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))
# # Convert to sf object
# pop_by_catchments_sf <- sf::st_as_sf(pop_by_catchments)
return(out = pop_by_catchments)
}
# Invoking the function ---------------------------------------------------------------------------------------
malaria_pop_by_catchment_2017 <- extract.pop.values(kasungu_population_2017, malaria_by_catchment)
malaria_pop_by_catchment_2018 <- extract.pop.values(kasungu_population_2018, malaria_by_catchment)
malaria_pop_by_catchment_2019 <- extract.pop.values(kasungu_population_2019, malaria_by_catchment)
malaria_pop_by_catchment_2020 <- extract.pop.values(kasungu_population_2020, malaria_by_catchment)View population maps
Estimated total number of people within health facility catchment areas.
# Function to create maps of estimated population by catchment areas ------------------------------------------------
create.population.map <- function(catchment.area, variable = "pop", title, legend.title = "Estimated \n population"){
# estimated population map
# catchment.area: estimated population layer from nachulu function
# variable: variable name (as character, in qoutes)
# title: map title in quotes
# legend.title: legend title in qoutes
# returns:
# a tmap-element (plots a map)
tm_shape(catchment.area)+
tm_fill(col = variable,
breaks = c(0, 13000, 19000, 27000, 35000, 70000, 140000, 200000),
palette = "YlOrBr",
title = legend.title)+
tm_borders(col = "grey",
lwd = 0.4)+
tm_layout(legend.position = c("left","bottom"),
legend.text.size = 0.6,
legend.title.size = 0.8,
frame = FALSE)+
tm_credits(title,
position = c(0.1, 0.8),
size = 1)+
tm_layout(legend.outside = TRUE)
}
# Invoking the function -----------------------------------------------------------------------------------------
pop_by_catchment_2017 <- create.population.map(malaria_pop_by_catchment_2017, title = "2017")
pop_by_catchment_2018 <- create.population.map(malaria_pop_by_catchment_2018, title = "2018")
pop_by_catchment_2019 <- create.population.map(malaria_pop_by_catchment_2019, title = "2019")
pop_by_catchment_2020 <- create.population.map(malaria_pop_by_catchment_2020, title = "2020")
tmap::tmap_arrange(pop_by_catchment_2017, pop_by_catchment_2018,
pop_by_catchment_2019, pop_by_catchment_2020, ncol = 2)
Fig. 2: Estimated population in Kasungu health facility catchment areas
Calculate the expected number of cases for each catchment area
Expected number of dry season malaria cases are calculated under the assumption that there is no spatial variation in risk, that is, no difference in infection rates between the catchment areas.
# Calculate expected malaria cases --------------------------------------------------------------
expected_malaria_2017 <- malaria_pop_by_catchment_2017 %>%
dplyr::rename(observed_2017 = dr_2017,
pop_2017 = pop) %>%
dplyr::mutate(expected_2017 = round(sum(observed_2017)/sum(pop_2017, na.rm = TRUE)*pop_2017))
expected_malaria_2018 <- malaria_pop_by_catchment_2018 %>%
dplyr::rename(observed_2018 = dr_2018,
pop_2018 = pop) %>%
dplyr::mutate(expected_2018 = round(sum(observed_2018)/sum(pop_2018, na.rm = TRUE)*pop_2018))
expected_malaria_2019 <- malaria_pop_by_catchment_2019 %>%
dplyr::rename(observed_2019 = dr_2019,
pop_2019 = pop) %>%
dplyr::mutate(expected_2019 = round(sum(observed_2019)/sum(pop_2019, na.rm = TRUE)*pop_2019))
expected_malaria_2020 <- malaria_pop_by_catchment_2020 %>%
dplyr::rename(observed_2020 = dr_2020,
pop_2020 = pop) %>%
dplyr::mutate(expected_2020 = round(sum(observed_2020)/sum(pop_2020, na.rm = TRUE)*pop_2020))
# Calculate Standard Morbidity Ratio (SMR) ---------------------------------------------------------
SMR_2017 <- expected_malaria_2017 %>%
dplyr::mutate(SMR = round(observed_2017/expected_2017,2)) %>%
dplyr::select(Names, pop_2017, observed_2017, expected_2017, SMR) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))
# SMR_2017$SMR[which(SMR_2017$Names == "Newa Mpasazi Health Centre")] <- 0
SMR_2018 <- expected_malaria_2018 %>%
dplyr::mutate(SMR = round(observed_2018/expected_2018, 2)) %>%
dplyr::select(Names, pop_2018, observed_2018, expected_2018, SMR) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))
SMR_2019 <- expected_malaria_2019 %>%
dplyr::mutate(SMR = round(observed_2019/expected_2019, 2)) %>%
dplyr::select(Names, pop_2019, observed_2019, expected_2019, SMR) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))
SMR_2020 <- expected_malaria_2020 %>%
dplyr::mutate(SMR = round(observed_2020/expected_2020, 2)) %>%
dplyr::select(Names, pop_2020, observed_2020, expected_2020, SMR) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))
SMR_table_2017 <- SMR_2017 %>%
kable %>%
kableExtra::kable_styling(full_width = FALSE)
SMR_table_2018 <- SMR_2018 %>%
kable %>%
kableExtra::kable_styling(full_width = FALSE)View observed and expected dry season malaria cases
# Function to create maps of observed and expected dry season malaria cases
create.malaria.map <- function(malaria.data, variable = NA, title = NA, legend.title = NA){
# observed and expected malaria incidence map
# malaria.data: data frame containing observed and expected malaria cases
# variable: variable name (as character, in qoutes e.g. "observed")
# title: map title in quotes
# legend.title: legend title in qoutes
# returns:
# a tmap-element (plots a map)
tm_shape(malaria.data)+
tm_fill(col = variable,
breaks = c(0, 500, 1000, 2500, 5000, 10000, 15000, 25000, 35000),
palette = "YlOrRd",
title = legend.title)+
tm_borders(lw = 0.3)+
tm_layout(legend.position = c(0.1,"bottom"),
legend.text.size = 0.5,
legend.title.size = 0.7,
frame = FALSE)+
tm_credits(title,
position = c(0.1, 0.8),
size = 1)+
tm_layout(legend.outside = TRUE)
}
# Invoking the function
# 2017 observed and expected malaria cases -------------------------------------------------------
observed_malaria_2017_map <- create.malaria.map(malaria_pop_by_catchment_2017,
variable = "dr_2017",
title = "2017",
legend.title = "Observed malaria")
expected_malaria_2017_map <- create.malaria.map(expected_malaria_2017,
variable = "expected_2017",
title = "2017",
legend.title = "Expected malaria")
# 2018 observed and expected malaria cases -------------------------------------------------------
observed_malaria_2018_map <- create.malaria.map(malaria_pop_by_catchment_2018,
variable = "dr_2018",
title = "2018",
legend.title = "Observed malaria")
expected_malaria_2018_map <- create.malaria.map(expected_malaria_2018,
variable = "expected_2018",
title = "2018",
legend.title = "Expected malaria")
# 2019 observed and expected malaria cases ---------------------------------------------------------
observed_malaria_2019_map <- create.malaria.map(malaria_pop_by_catchment_2019,
variable = "dr_2019",
title = "2019",
legend.title = "Observed malaria")
expected_malaria_2019_map <- create.malaria.map(expected_malaria_2019,
variable = "expected_2019",
title = "2019",
legend.title = "Expected malaria")
# 2020 observed and expected malaria cases ---------------------------------------------------------
observed_malaria_2020_map <- create.malaria.map(malaria_pop_by_catchment_2020,
variable = "dr_2020",
title = "2020",
legend.title = "Observed malaria")
expected_malaria_2020_map <- create.malaria.map(expected_malaria_2020,
variable = "expected_2020",
title = "2020",
legend.title = "Expected malaria")
# Layout maps
tmap::tmap_arrange(observed_malaria_2017_map, expected_malaria_2017_map,
observed_malaria_2018_map, expected_malaria_2018_map,
observed_malaria_2019_map, expected_malaria_2019_map,
observed_malaria_2020_map, expected_malaria_2020_map, ncol = 2)
Fig 3: Observed and expected malaria incidence by health facility catchment area, Kasungu
View SMR by catchment
# Function to create maps of Standard Morbidity Rate by catchment --------------------------
create.smr.map <- function(smr.data, variable = "SMR", title = NA, legend.title = "SMR"){
# SMR by catchment map
# smr.data: sf polygon object containing SMR by catchment
# variable: variable name (as character, in qoutes)
# title: map title in quotes
# legend.title: legend title in qoutes
# returns:
# a tmap-element (plots a map)
tm_shape(smr.data)+
tm_fill(col = variable,
breaks = c(0, 0.5, 1, 1.5, 5, 10, 50, 100, 155, 200),
palette = "-magma",
title = legend.title)+
tm_borders(lw = 0.3)+
tm_layout(legend.position = c(0.1,"bottom"),
legend.text.size = 0.5,
legend.title.size = 0.7,
frame = FALSE)+
tm_credits(title,
position = c(0.1, 0.8),
size = 1)+
tm_layout(legend.outside = TRUE)
}
# Invoking function -------------------------------------------------------------------------
SMR_2017_map <- create.smr.map(SMR_2017, title = "2017")
SMR_2018_map <- create.smr.map(SMR_2018, title = "2018")
SMR_2019_map <- create.smr.map(SMR_2019, title = "2019")
SMR_2020_map <- create.smr.map(SMR_2020, title = "2020")
# Layout maps
tmap::tmap_arrange(SMR_2017_map, SMR_2018_map, SMR_2019_map, SMR_2020_map, ncol = 2)
Fig. 4: Standard Morbidity Rate by health facility catchment